Computing the current density j in the coil¶
\(\def\curl{\operatorname{curl}}\def\Curl{\operatorname{Curl}}\def\div{\operatorname{div}}\) Before we solve the non-linear magnetostatic problem, we first need to find the current \(j\) flowing inside the coil \(\Omega_c\).
The first property is that \(j\) is solenoidal, i.e. \(\div(j)=0\). Further, we presume the existence of a potential \(\phi\) with \(j=-\sigma\nabla\phi\), where \(\sigma\) denotes the connectivity. As for the boundary conditions, we prescribe \(n\cdot j = -\sigma\partial_n\phi=0\) on the exterior boundary \(\Gamma_{ex}\), and \(\phi=1\) and \(\phi=0\) on the inflow \(\Gamma_{in}\) and outflow boundary \(\Gamma_{out}\), respectively. Altogether, we have
However, in our case, the coil is a loop with no inflow or outflow boundary. This is where the face “coil_cut_1” we defined in the geometry comes into play. For this purpose, we have to introduce “fictitious” points and introduce a clone of the face “coil_cut_1” in order to be able to prescribe the necessary boundary conditions.
Ok so how do we even solve this? Weak formulation leads to
The boundary condition \(\partial_n\phi=0\) on \(\Gamma_{ex}\) is natural, while the Dirichlet conditions are essential and have to included in the space. Further, since we solve the problem only on the coil, the other DOFs have to be eliminated from the final system.
We begin by loading the geometry and the mesh generated in the previous document
[49]:
%%capture
%run TEAM_13_geometry.ipynb
Create the MESH object from the mesh generated by netgen
[50]:
import sys
sys.path.insert(0,'../../../') # adds parent directory
import pde
MESH = pde.mesh3.netgen(geoOCCmesh)
MESH
[50]:
np:4132, nt:23369, nf:3059, ne:720, nf_all:46966, ne_all:27728
The piece of code below duplicates the face, as described above, and generates a new MESH object.
[51]:
import numpy as np
face_index = pde.tools.getIndices(MESH.regions_2d,'coil_cut_1')
faces = MESH.f[MESH.BoundaryFaces_Region == face_index,:3]
new_faces = faces.copy()
points_to_duplicate = np.unique(faces.ravel())
new_points = np.arange(MESH.np, MESH.np+points_to_duplicate.size)
actual_points = MESH.p[points_to_duplicate,:]
t_new = MESH.t[:,:4].copy()
f_new = MESH.f[:,:3].copy()
p_new = MESH.p.copy()
for i,pnt in enumerate(points_to_duplicate):
# append point to list
p_new = np.vstack([p_new,p_new[pnt,:]])
# finding tets coordinates containing the ith point to duplicate
tets_containing_points = np.argwhere(t_new[:,:4]==pnt)[:,0]
for _,j in enumerate(tets_containing_points):
#check if tet is left
if MESH.mp_tet[j,0]<0:
t_new[j,t_new[j,:]==pnt] = MESH.np + i
# finding faces containing the points
faces_containing_points = np.argwhere(f_new[:,:3]==pnt)[:,0]
for _,j in enumerate(faces_containing_points):
#check if face is left
if 1/3*(p_new[f_new[j,0],0] + p_new[f_new[j,1],0] + p_new[f_new[j,2],0])<0:
f_new[j,f_new[j,:]==pnt] = MESH.np + i
# print(faces_containing_points)
t_new = np.c_[t_new,MESH.t[:,4]]
for i,j in enumerate(faces.ravel()):
new_faces.ravel()[i] = new_points[points_to_duplicate==j][0]
new_faces = np.c_[new_faces,np.tile(MESH.f[:,3].max()+1,(new_faces.shape[0],1))]
f_new = (np.r_[np.c_[f_new,MESH.f[:,3]],new_faces]).astype(int)
regions_2d_new = MESH.regions_2d.copy()
regions_2d_new.append('new')
identifications = (np.c_[points_to_duplicate,new_points]).astype(int)
# stop
MESH = pde.mesh3(p_new,MESH.e,f_new,t_new,MESH.regions_3d,regions_2d_new,MESH.regions_1d,identifications = identifications)
MESH
[51]:
np:4149, nt:23369, nf:3076, ne:720, nf_all:47028, ne_all:27806
As we can see, additional points and faces have been created to accommodate the additional interface.
Next, we proceed to solve the problem in \(H^1\). In weak form, we have:
Find \(\phi\in V^* = \{u\in H^1\;:\; u|_{\Gamma_{out}}=0 \text{ and } u|_{\Gamma_{in}}=1\}\) such that \((\sigma\nabla\phi,\nabla v) = 0\), for all \(v\in H^1\)
After homogenization, which involves splitting the solution \(\phi = \phi_* + \phi_0\), where \(\phi_*\in V^*\), we can instead solve:
Find \(\phi_0\in H^1\setminus\Gamma_{in,out}\) such that \((\sigma\nabla\phi_0,\nabla v) = -(\sigma\nabla\phi_*,\nabla v)\), for all \(v\in H^1\setminus\Gamma_{in,out}\)
After using conforming \(P_1\) continuous finite elements, we solve the discretized system in the following way: we split the discrete solution \(\phi = R^T_{in}\phi_{in} + R^T_{out}\phi_{out} + R^T_{int}\phi_{int} = R^T_{in}\phi_{in} + R^T_{int}\phi_{int}\), where \(R^T_*\) matrices which assign the correct indices to the boundaries and the interior degrees of freedom. Assume the stiffness matrix is denoted by \(K\) and the right-hand side is \(r\). Then, purely in \(H^1\), the system is of the form \(K\phi=f\). Plugging in the splitting for \(\phi\) leads to the system \(KR_{int}^T\phi_{int} = -KR_{in}^T\phi_{in}\). By multiplying from the left with \(R_{int}\), we obtain the square system
[52]:
order = 1
D = pde.int.assemble3(MESH, order = order)
DB = pde.int.assembleB3(MESH, order = order)
N1,N2,N3 = pde.int.assembleN3(MESH, order = order)
unit_coil = pde.int.evaluate3(MESH, order = order, coeff = lambda x,y,z : 1+0*x, regions = 'coil')
###########################################################################
phi_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'M', order = order)
dphix_H1, dphiy_H1, dphiz_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = order)
phiB_H1 = pde.h1.assembleB3(MESH, space = 'P1', matrix = 'M', shape = phi_H1.shape, order = order)
R0, RS0 = pde.h1.assembleR3(MESH, space = 'P1', faces = 'new,coil_cut_1')
R1, RS1 = pde.h1.assembleR3(MESH, space = 'P1', faces = 'coil_cut_1')
M = phi_H1 @ D @ unit_coil @ phi_H1.T
K = dphix_H1 @ D @ unit_coil @ dphix_H1.T +\
dphiy_H1 @ D @ unit_coil @ dphiy_H1.T +\
dphiz_H1 @ D @ unit_coil @ dphiz_H1.T
r = -RS0 @ K @ R1.T @ (1+np.zeros(R1.shape[0]))
K = RS0 @ K @ RS0.T
RZ = pde.tools.removeZeros(K)
K = RZ @ K @ RZ.T
r = RZ @ r
sigma = 1#58.7e6
from sksparse.cholmod import cholesky as chol
p = chol(sigma*K).solve_A(r)
potential_H1 = RS0.T @ RZ.T @ p + R1.T @ (1+np.zeros(R1.shape[0]))
jx_L2 = -(dphix_H1.T@potential_H1)*unit_coil.diagonal()
jy_L2 = -(dphiy_H1.T@potential_H1)*unit_coil.diagonal()
jz_L2 = -(dphiz_H1.T@potential_H1)*unit_coil.diagonal()
dphix_H1_P0, dphiy_H1_P0, dphiz_H1_P0 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = 0)
unit_coil_P0 = pde.int.evaluate3(MESH, order = 0, coeff = lambda x,y,z : 1+0*x, regions = 'coil')
jx_L2_P0 = -(dphix_H1_P0.T@potential_H1)*unit_coil_P0.diagonal()
jy_L2_P0 = -(dphiy_H1_P0.T@potential_H1)*unit_coil_P0.diagonal()
jz_L2_P0 = -(dphiz_H1_P0.T@potential_H1)*unit_coil_P0.diagonal()
Instead of solving in \(H^1\), we can instead solve the mixed problem directly by employing \(H(\div)\) conforming finite elements.
[53]:
order = 1
phix_Hdiv, phiy_Hdiv, phiz_Hdiv = pde.hdiv.assemble3(MESH, space = 'RT0', matrix = 'M', order = order)
divphi_Hdiv = pde.hdiv.assemble3(MESH, space = 'RT0', matrix = 'K', order = order)
phi_L2 = pde.l2.assemble3(MESH, space = 'P0', matrix = 'M', order = order)
D = pde.int.assemble3(MESH, order = order)
M_Hdiv_coil_full = phix_Hdiv @ D @ unit_coil @ phix_Hdiv.T +\
phiy_Hdiv @ D @ unit_coil @ phiy_Hdiv.T +\
phiz_Hdiv @ D @ unit_coil @ phiz_Hdiv.T
C_Hdiv_L2 = divphi_Hdiv @ D @ unit_coil @ phi_L2.T
R1, RS1 = pde.hdiv.assembleR3(MESH, space = 'RT0', faces = 'coil_face')
M_Hdiv_coil_full = RS1 @ M_Hdiv_coil_full @RS1.T
C_Hdiv_L2 = RS1 @ C_Hdiv_L2
import scipy.sparse as sp
AA = sp.bmat([[M_Hdiv_coil_full, C_Hdiv_L2],
[C_Hdiv_L2.T, None]])
RZdiv = pde.tools.removeZeros(AA)
AA = RZdiv @ AA @ RZdiv.T
phiB_Hdiv = pde.hdiv.assembleB3(MESH, space = 'RT0', matrix = 'phi', shape = phix_Hdiv.shape, order = order)
unit_coil_B = pde.int.evaluateB3(MESH, order = order, coeff = lambda x,y,z : 1+0*x, faces = 'new').diagonal()
DB = pde.int.assembleB3(MESH, order = order)
rhs = unit_coil_B @ DB @ phiB_Hdiv.T
rhs = np.r_[RS1@rhs,np.zeros(MESH.nt)]
rhs = RZdiv @ rhs
xx = sp.linalg.spsolve(AA,rhs)
potential_L2 = (RZdiv.T@xx)[-MESH.nt:]
j_hdiv = RS1.T@(RZdiv.T@xx)[:-MESH.nt]
jx_hdiv = (phix_Hdiv.T@j_hdiv)*unit_coil.diagonal()
jy_hdiv = (phiy_Hdiv.T@j_hdiv)*unit_coil.diagonal()
jz_hdiv = (phiz_Hdiv.T@j_hdiv)*unit_coil.diagonal()
##############################################################################
unit_coil_P0 = pde.int.evaluate3(MESH, order = 0, coeff = lambda x,y,z : 1+0*x, regions = 'coil')
phix_Hdiv_P0, phiy_Hdiv_P0, phiz_Hdiv_P0 = pde.hdiv.assemble3(MESH, space = 'RT0', matrix = 'M', order = 0)
jx_hdiv_P0 = (phix_Hdiv_P0.T@j_hdiv)*unit_coil_P0.diagonal()
jy_hdiv_P0 = (phiy_Hdiv_P0.T@j_hdiv)*unit_coil_P0.diagonal()
jz_hdiv_P0 = (phiz_Hdiv_P0.T@j_hdiv)*unit_coil_P0.diagonal()
##############################################################################
[54]:
grid = pde.tools.vtklib.createVTK(MESH)
pde.tools.vtklib.add_H1_Scalar(grid, potential_H1, 'potential_H1')
pde.tools.vtklib.add_L2_Scalar(grid, potential_L2, 'potential_L2')
pde.tools.vtklib.add_L2_Vector(grid,jx_L2_P0,jy_L2_P0,jz_L2_P0,'J_L2')
pde.tools.vtklib.add_L2_Vector(grid,jx_hdiv_P0,jy_hdiv_P0,jz_hdiv_P0,'J_HDIV')
pde.tools.vtklib.writeVTK(grid, 'current_density.vtu')
[55]:
import pyvista as pv
mesh = pv.read('current_density.vtu')
mesh
[55]:
| Header | Data Arrays | ||||||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
[56]:
mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,1])
p = pv.Plotter()
threshed.set_active_scalars("potential_H1")
p.add_mesh(threshed, style='surface', opacity = 0.4, label=None)
mesh.set_active_vectors("J_L2")
arrows = mesh.glyph(scale="J_L2", orient=True, tolerance=0.03, factor=9500.0)
p.add_mesh(arrows, color="black")
p.camera_position = [(0, 0, 600),(0, 0, 0),(0, 0, 0)]
p.show(jupyter_backend='html')
[57]:
mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,1])
p = pv.Plotter()
threshed.set_active_scalars("potential_L2")
p.add_mesh(threshed, style='surface', opacity = 0.4, label=None)
mesh.set_active_vectors("J_HDIV")
arrows = mesh.glyph(scale="J_HDIV", orient=True, tolerance=0.03, factor=9500.0)
p.add_mesh(arrows, color="black")
p.camera_position = [(0, 0, 600),(0, 0, 0),(0, 0, 0)]
p.show(jupyter_backend='html')